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ABSTRACT 

We investigate the coupled formation and evolution of galaxies and their embedded supermassive black holes 
using state-of-the-art hydrodynamic simulations of cosmological structure formation. For the first time, we self- 
consistently follow the dark matter dynamics, radiative gas cooling, star formation, as well as black hole growth 
and associated energy feedback processes, starting directly from initial conditions appropriate for the ACDM 
cosmology. Our modeling of the black hole physics is based on an approach we have recently developed and 
tested in simulations of isolated galaxy mergers. Here we apply the same model in cosmological simulations to 
examine: (i) the predicted global history of black hole mass assembly in galaxies, (ii) the evolution of the local 
black hole-host mass correlations and (iii) the conditions that allow rapid growth of the first quasars, as well 
as the properties of their hosts and descendants today. We find that our simulations produce a total black hole 
mass density pbh — 2 x 10^ Mq Mpc^"' by z = 0, in good agreement with observational estimates. The black 
hole accretion rate density, p^^, peaks at lower redshift and evolves more strongly at high redshift than the star 
formation rate density, p^, with an approximate scaling as Pbh/p* oc (1 + z)^* at z > 3. On the other hand, 
the ratio Pbh/p* ^ (1 +z)^"^ of black hole to stellar mass densities shows only a moderate evolution at low 
redshifts z < 3. For the population of galaxies identified in the simulations at z = 1 we find strong correlations 
between black hole mass and velocity dispersion or mass of the stellar systems. The predicted correlations 
agree well with the measured local Mbh — cr and Mbh — A^* relationships, but also suggest a weak evolution 
with redshift in the normalization, and in particular the slope. However, the magnitude of this effect is sensitive 
to the range of masses being probed. For stellar masses of > 3 x 10'*\ we predict a trend of increasing 
M^/Mbh with redshift, in agreement with recent direct estimates of the BH to host stellar mass ratio at high 
redshift and the conjecture that a more fundamental relation (a BH fundamental plane) should involve both 
and cr. We find that our simulation models can also produce quite massive black holes at high redshift, 
as a result of extended periods of exponential growth in relatively isolated, rare regions that collapse early 
and exhibit strong gas inflows. Interestingly, when followed to their descendants, these first supermassive BH 
systems are not necessarily the most massive ones today, since they are often overtaken in growth by quasars 
that form later. 

Subject headings: quasars: general — galaxies: formation — galaxies: active — galaxies: evolution — cos- 
mology: theory — hydrodynamics 



1. INTRODUCTION 

Following the discovery of quasars dSchmidtl 119631: 
iGreenstein & Matthewsl [l963l) it was suggested that super- 
massive black holes (10^ — lO' M0) lie at the centers of 
galaxies, and that the quasar activity is fueled by the re- 
lease of gravitational energy from their accreted matter The 
remnants of quasar phases at early times are probably the 
supermassive black holes found at the centers of galaxies 
in our local Universe. Interestingly, the properties of these 
supermassi ve black ho les are tightly coupled to the mass 
jMagorrian et al.lflQQSl) and velocity dispersion of their host 
galaxies, as manifested in the Mbh ~ rela tion of spheroids 
dFen-arese & Merritt 2000: Gebhardt et all I2OOO1) . In addi- 
tion, the black hole mass is correlated with the concentra- 
tion or Sers ic index dGraham & Driveij|2006l) . Most recently, 
iHopkins eT al. (2007al) have shown that these various correla- 
tions are not independent, and can be understood as projec- 
tions of a "black hole fundamental plane" (BHFP), similar to 
that describing properties of elliptical galaxies. 
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The existence of highly luminous quasars also constrains 
the formation and evolution of massive galaxies and the epoch 
of reionization. Quasars with inferred black hole masses in 
excess of 10^ have now been discovered out to z 6 
dFan et al.ll2003h . indicating an early formation time for black 
holes and galaxy spheroids and posing a significant challenge 
for theoretical models of high-redshift quasar and galaxy for- 
mation. 

This growing observational evidence, drawn from local 
galaxies to high redshift quasars, argues for a close connec- 
tion between the formation and evolution of galaxies and of 
their central supermassive black holes. However, the physical 
nature of this relationship has yet to be understood in detail. 
Indeed, there are significant gaps in our observational and the- 
oretical knowledge of the history of black hole formation and 
evolution in galaxies. 

For example, current velocity dispersion measurements 
are inconclusive about the important question whether the 
tight scaling relations evolve with redshift (Woo et al. 2006t 
Shields et al. 2006), or are essentially invariant (Shields et al] 
2003ir as a function of time. We note that some evolution 
in the ratio of black h ole to halo mass i s suggest ed by clus- 
tering constraints (e.g lAdelberger & Stei del 2005: lLidz et al.l 
12006), but these measurements do not directly address the 
relationship between black hole and properties of the lumi- 
nous host galaxy. More relevant are comparisons of the black 
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hole mass inferred from quasar observations to the host stel- 
lar m ass, both obseryationa lly (lMerlonill2004l) and theoreti- 
cally dHopkins et al.l 12006 a'). which indicate an evolution in 
e.g. the Magorrian relation, in the sense that black holes are 
more massive relative to luminous spheroids at high redshifts 
than at z 0. 

Theoretical studies of the co-evolution of black holes 
and galaxies have so far mostly used so-called semi- 
analytical modeling (e.g. IjC auffmann & Haehnelt 2000; 
Cattaneo et al. 1999 1 IWyithe & Loebl l2003t iVolonteri et alj 
20031; iDi Matteo etal.1 120031: |i 
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2006t " lDe Lucia et al.l 12006^ .Malbon et al.1 l2007h of galaxy 
formation, in which the growth of galaxies and their 



embedded black holes is followed with simple physical pa- 
rameterizations on top of dark matter merging history trees. 
Many of these models assume that quasar activity is triggered 
by major galaxy mergers, motivated by hydrodynamical sim- 
ulations that have shown that gravitational tidal fields during 
major mergers of gas rich galaxies produce strong gas inflows 
(jBarnes & Hernquist 1991, 1996), which le ad to a burst of 
nuclear star formation dMihos & Hernquist i 19961) and are 
likely the prerequisite for rapid black hole growth and quasar 
activity. Nearby quasars are indeed pref erentially found in 
tidally disturbed objects (e.g Jogee' 2004), corroborating the 
importance of galaxy interactions and mergers for major 
black hole growth. 

Many theoretical explanations for the observed correlations 
between galaxy properties and black hole mass rely on some 
form of self-regulated growth of the BHs. For example, it has 
been suggested that the central black holes grow until they 
release sufficient energ y to unbind the gas that feeds them 
from their host gal a xy jCiotti & Ostrikerl l 19971; ISilk & ReesI 
[19981 iFabianI [19991: iWvithe & Loeb 20Q1. We recently ex- 
plored such a local energy feedback for the first time with self- 
consistent, fully three-dimensional hydrodynamic simulations 
of galaxies (IDi Matteo et al.ll2005):[Springel et al.ll2005bl) that 
include a treatment of accretion on supermassive black holes 
and their associated energy feedback. These simulations have 
demonstrated that the fundamental BH-host correlation in- 
cluding the Mbh — cr relation can indeed b e reproduced in 
feedback-regulated models of BH growth dDi Matteo et alj 
[2005.) . in accordance with theoretical conjectures. At the 
same time, the dynamical coupling in the simulations of hy- 
drodynamical gas inflow, star formation, black hole growth 
and associated feedback processes gives them substantial pre- 
dictive power well beyond that of simplified analytical and 
semi-analytical model s . Besides the Mbh ~ or Mbh — A^* 
dDi Matteo et al.ll2005l: [Robertson et al.ll2006bh relationships, 
the simulation models can for example predict the detailed 
properties of the spheroidal galaxies forming in major merg- 
ers and how they correlate with the BH masses. In fact, 
they suggest the existenc e of a fundamental p lane relation 
for BHs (Mbh oc a^ '^R^/ dHopkins et al.ll2007al) . provide an 
explanation for the red colors of massive elliptical galaxies 
((Springel et al. 2005 ^, and describe the fundamental plane of 
elliptical galaxies ( Robertson et al.l l2006a). They also suggest 
luminosity-dependent quasar lifetimes, leading to a new in- 
terpretation for the origin of the quasar luminosity function 
and its evolution over cosmic history dHopkins et al.l [20051 
[l306a). 

In the present paper, we extend these earlier studies by car- 
rying out fully cosmological hydrodynamic simulations of 
the ACDM model that jointly follow the growth of galax- 



ies and supermassive black holes, as well as their associ- 
ated feedback processes. Our approach is based on the same 
methodology that we have developed and applied in the high- 
resolution simulations of galaxy mergers, augmented with a 
suitable mechanism to seed emerging new dark matter halos 
with a small black hole that can then grow by gas accretion 
later on. While much more restricted in numerical resolution 
than simulations of individual galaxy mergers, our modeling 
of star formation and black hole physics in terms of a sub- 
resolution treatment provides quite accurate results already at 
comparatively coarse resolution, an important prerequisite for 
attempting to model these processes in cosmological simu- 
lations. Nevertheless, numerical resolution is clearly an im- 
portant limitation of our cosmological results, an aspect that 
we will discuss in more detail where appropriate. With this 
caveat in mind, we would like to stress however that the un- 
ambiguous initial conditions of direct cosmological simula- 
tion make them in principle the most powerful and accurate 
tool for studying the interplay of galaxy formation and black 
hole growth. Our aim is therefore to examine how well our 
current model for treating BH physics in simulations does in 
present state-of-the-art hydrodynamical calculations of cos- 
mic structure formation, and what we can learn from them 
to advance our theoretical understanding of the co-evolution 
of galaxies and supermassive black holes. In this study we 
shall focus on basic properties of the black hole population, 
like the evolution of the cosmic BH mass density, and the cor- 
relations between B H masses and host galaxy properties. In 
ISiiacki et al.l d2007l) we also study an extension of our feed- 
back model with a 'radio mode' that is active at low accretion 
rates and is distinct from the normal quasar activity, allowing 
us to study the formation of AGN in rich galaxy clusters at 
low redshifts, where it is likely to be important. 

This paper is structured as follows. In §2 we describe our 
simulation set and the numerical modeling adopted for the 
ISM, star formation and gas accretion onto black holes. In 
§3 we present our results for the evolution of the global black 
hole mass density and compare it to the cosmic history of star 
formation and the evolution of the stellar density. In §4 we 
examine our results for the fundamental BH/host correlations 
measured from the simulation and for their evolution from 
high to low redshift. Finally, we summarize and discuss our 
findings in §5. 

2. METHODOLOGY 

2.1. Numerical code 

In this study we focus on a ACDM cosmological model 
with parameters chosen according to the first year results 
from the WiUdnson Microwave Anisotropy Probe (WMAPl; 
ISpergel et aLll2003h . f^o = 0.3, = 0.7, Hubble constant 
7^0 = 100/!kms-i Mpc"' with h = 0.7 and a scale invari- 
ant primordial power spectrum with index « = 1, with a 
normalization of the amplitude of fluctuations erg = 0.9. 
We use a significantly extended version of the parallel cos- 
mological TreePM-SPH code GADGET2 (Springel 2005) to 
evolve a realization of ACDM initial conditions from high 
to low redshift. The combination of a high-resolution grav- 
itational solver with individual and adaptive timesteps allows 

The largest simulation presented here had already been started by the 
time the updated third year constraints have become available (WMAP3; 
Spergel et al. 2006). We comment on effects on the growth of the halo mass 
fu nction owing to t he lower amplitude o f fluctuations, erg implied by WMAP3 
in ILi et all <2007h : [Sliacld et alj flOOTi) 
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this code to bridge a large dynamic range both in length- and 
timescales. Gas dynamics is followed with the Lagrangian 
smoo thed particle hydrodynamics (SPH) (e.g iMonaghanI 
II992I) technique, which we employ in a formulation that 
manifestly conserves energy and entropy, despit e the use of 
fully adaptive SPH smoothing lengths ( S pringel & HernquistI 
|2002|). Radiative coolin g and heating processes are computed 
as in iKatz et alj (Il996h . with a spatially uniform photoioniz- 
ing UV background that is imposed externally. 

Within cosmological (or galaxy-sized) numerical simula- 
tions, it is presently (and for some time to come) not feasible 
to follow the physics of star formation and black hole accre- 
tion from first principles down to scales of individual stars or 
black holes. Any numerical model of galaxy formation there- 
fore needs to make substantial approximations for some of the 
relevant physics on unresolved scales. 

For modeling star formation and its associated supernova 
feedback we use the sub-resolution multiphase model for 
the int erstellar medium developed by ISpringel & HernquistI 
((2003 a). In this model, a thermal instability is assumed to 
operate above a critical density threshold pth, producing a two 
phase medium consisting of cold clouds embedded in a ten- 
uous gas at pressure equilibrium. Stars form from the cold 
clouds, and short-lived stars supply an energy of 10^' ergs to 
the surrounding gas as supernovae. This energy heats the dif- 
fuse phase of the ISM and evaporates cold clouds, thereby 
establishing a self-regulation cycle for star formation, pth is 
determined self-consistently in the model by requiring that 
the equation of state (EOS) is continuous at the onset of 
star formation. The cloud evaporation process and the cool- 
ing function of the gas then determine the temperatures and 
the mass fractions of the two hot and cold phases of the 
ISM, such that the EOS of the model can be directly com- 
puted as a function of density. The latter is encapsulating 
the self-regulated nature of star formation owing to super- 
novae feedback in a s imple model for a multiphase ISM. As 
in the Springel & HernquistI (l2003a) model we have included 
a model for supernova-driven galactic winds with an initial 
wind speed of v ^ 480 km s^' . 

For the parameter settings adopted here, the model repro- 
duces the observed star formation rate surface densities in iso- 
lated spiral galaxies (Kennicutti ll989i Il998l) . Using a large 
number of nested cosmological simulations, the approach we 
adopt (and the parameters we use) has also been shown to 
lead to a numerically converged estimate for the cosmic star 
formation history of the universe that a grees reasonably well 
with low redshif t observations (Springe l & Hernq uist 2003bt 
iHernquist & Springel 200 3). For the modeUng of BH accre- 
tion and feedback we adopt a similar strategy as for star for- 
mation, which we discuss next. 

2.2. Accretion and Feedback from Supermassive black holes 

As for star formation, we adopt a sub-resolution model to 
capture the main features of accretion and assoc iated feedback 
on supermassive black holes (as introduced in ISpringel et al.l 
l2005bl:lDrMatteo et al.l2005l) . To this end, we represent black 
holes by collisionless 'sink' particles that can grow in mass by 
accreting gas from their immediate environments, or by merg- 
ing with other black holes. We estimate the gas accretion rate 
onto a black hole using a Bondi-Hoyle-Lyttleton parameteri- 
zation (Bondi 1952; Bond i & Hoylelll944l:lHoyle & Lyttletoiil 
Il939l) . In this description, the accretion rate onto the black 
hole is given by Mb = 47r [G^M|jj/9]/(c^ + v^)-'/^ where p, 
C.5 are the density and sound speed of the ISM gas, respec- 



tively, and V is the velocity of the black hole relative to the 
surrounding gas. In this study we also allow the black hole to 
accrete at mildly super-Eddington values and impose a maxi- 
mum allowed accretion rate equal to 2x Eddington rate, MEdd- 
As we shall see such conditions are generally achieved at high 
redshift. We note that the detailed relativistic accretion flow 
onto the black hole is unresolved in our simulations, but if 
the limitating factor for rapid growth of BHs lies in the larger- 
scale gas distribution around the black hole, which is resolved, 
then the Bondi prescription should capture the dependence of 
the mean accretion rate onto the conditions of the gas in the 
region around the black hole. 

The radiated luminosity, L, is related to the gas accretion 
rate, Mbh, by the radiative efficiency 77 — Li-/ {Mbh c^) , which 
simply gives the mass to energy conversion efficiency, set by 
the amount of energy that can be extracted from the innermost 
stable orbit of an accretion disk around a black hole. We will 
adopt a fixed value o f er = 0.1, which is the mean value for a 
radiativelv efficient ( Shakura & Sunvaev 1973) accretion disk 
onto a Schwarzschild black hole. 

We assume that some small fraction Cf of the radiated lu- 
minosity Lr can couple to the surrounding gas in the form of 
feedback energy, viz. fifeed = CfLi In accordance with our pre- 
vious studies of galaxy merger simulations we take ef — 5%. 
This value governs the normalization of the Mbh — cr relation, 
and ef = 5% brings it in to agreement with current observations 
dPi Matteo et alfcoOSl) . We note that ef is effectively the only 
free parameter in our black hole model. After fixing it to re- 
produce the normalization of the observed Mbh ~ o" from the 
galaxy models we do not vary it in any of our cosmological 
simulations. 

For simplicity, we deposit the feedback energy isotropically 
in the region around the black hole. Lack of spatial resolution 
precludes us from considering mechanical modes of releas- 
ing the energy, e.g. in the form of a jet. However, we note 
that is is plausible that such other forms of energy are ther- 
malized eventually as well, and that the final impact of the 
feedback depends primarily on the total amount of energy re- 
leased and less on the form it is released in. This is likely 
to be true generally, provided that the energy (or momentum) 
is imparted to the surrounding gas on length scales small and 
time scales short compared with those that characterize the 
host galaxy. In that event, the impact of black hole feedback 
will be explosive in nature and, indeed, the blowout phase of 
evolution in our simulations is well-d escribed by a general- 
ized Sedov-Taylor blast- wav e solution (iHopkins et al.ll2006bl : 
iHopkins & Hernquistll2006l) . In any case, we emphasize that 
despite an isotropic release of the energy, the response of the 
gas can still be decidedly anisotropic, e.g. when a dense gas 
disk is present that channels the gas response into a collimated 
outflow. 

The idea that we follow with our feedback modeling here 
is the rapid accretion phases of BHs at times close to their 
critical growth phases. Such 'quasar' phases are typically rel- 
atively short-lived and require galaxy mergers to produce the 
strong gravitational tidal forcing necessary for sufficient nu- 
clear gas inflow rates. It is presently unclear whether the ac- 
cretion disks in such modes actually produce mechanical jets 
of release their feedback in another way. It is however evident 
that at low redshift (e.g.; z < 1 which we are not studying di- 
rectly here) that when BHs that are embedded in the hot gas 
atmospheres of large groups and clusters channel their energy 
in mechanical feedback in the form of relativistic jets, gen- 
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Fig. 1. — Projected baryonic density field in slices of thickness 5000/!"' kpc through our high resolution simulation, color-coded by temperature and with 
brightness proportional to the logaiithm of the gas density. Each panel shows the same region of space at different redshifts, as labeled. The circles mark the 
positions of the black holes, with a size that encodes the BH mass, as indicated in the top left panel. 
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TABLE 1 

Numerical parameters of cosmological simulations with BHs 



Run 


Boxsize 




N„ 


'«DM 


nigas 


e 


^end 












h-' 


Mq 


kpc 




D4 


33.75 


2 


x216^ 


2.75 X 10** 


4.24 


X 10' 


6.25 


0.00 


D6 (BHCosmo) 


33.75 


2 


X4863 


2.75 X 10' 


4.24 


X 10* 


2.73 


1.00 


E6 


50 


2 


X4863 


7.85 X lO' 


1.21 


X 10' 


4.12 


4.00 



erating buoyantly rising radio bubble. In lSiiacki et alj ( l2007h 
we consider this form of feedback mode. 

An important remaining question in our model concerns 
the ultimate origin of the BHs. Since our accretion rate es- 
timate can only grow a BH that already exists, we assume that 
a physical process that produces small seed BHs is operating 
sufficiently efficiently that effectively all halos above a certain 
threshold mass contain at least one such seed BH. Whether 
or not they can then grow to larger masses by gas accretion 
will be determined by the local gas conditions, as described 
in our model above. In order to achieve such a seeding at 
a technical level in cosmological runs we use an on-the-fly 
'friends-of-friends' group finder algorithm which is called at 
intervals equally spaced in the logarithm of the scale-factor 
a, with A logo — log 1.25. This provides the locations and 
mass of all halos in the simulation. If a halo is more massive 
than our threshold and does not contain any black hole yet, we 
endow it with one by converting its densest gas particle into 
a sink particle with a seed black hole mass of 10^/;^' M0. 
The further growth of the black hole sink can then proceed by 
gas accretion, at a rate that depends sensitively on the local 
conditions, or by mergers with other black hole sink parti- 
cles. The total cumulative black hole mass introduced in this 
way as seeds is negligible compared to the mass growth by 
gas accretion. We note that being able to run a fast, parallel 
'friends-of-friends' algorithm on the fly during simulations is 
an important technical prerequisite of our technique. 

Further motivation for this choice of seeding procedure is 
based on the currently proposed scenarios for the seed black 
holes in galaxies. To grow a supermassive black hole to a 
mass of ^ 10^ M0 in less than a billion years, as required 
by presence of the z = 6 SDSS quasars, may require (1) 
the catastrophic collapse of a supermassive star that forms a 
large i nitial black hole of mass 10"* - 10^ (ICair & ReesI 
Il984t iBromm & Loebl 120031: llegelman et al.1 120061) . or (2) 
alternatively, smaller black hole seeds (M ~ 10^ M©) may 
form from the first PopIII stars at 30 and grow exponen- 
tially from the n on tAbel et alJl2002HBromm & Larsonll2004l: 
lYoshidaetalJl2006l) . In our simulations, black hole seeds of 
mass M = 10^ Mq are introduced into galaxies as they 
initially reach Mhaio = 10'"/;^' Mq. This choice is a good 
approximation to what is expected for both of the hypotheses 
outlined above. For (1) this is roughly in the correct range; 
whereas for (2) Eddington growth predicts that the black hole 
has grown to roughly these values by the time of collapse of 
M 10"^/?^' Mq perturbations, which occurs at z 10 in our 
standard cold dark matter scenario. Additionally, although not 
required, this value of the initial black hole mass to galaxy ra- 
tio fits the observed relations at low redshift (iMagorrian et alJ 
Il998t iFerrarese & MerrittI 12000 ). It is important to note that 
the dominant growth of black holes always occurs in expo- 
nential Eddington phases induced by gas accretion so that our 
results are rather insensitive to the specific choice of the seed 
mass. 



When galaxies and their surrounding dark halos merge to 
form a single dark matter halo, their central black holes are 
also expected to merge eventually, so hierarchical black hole 
mergers contribute to the growth of the central black holes. 
Whether the forming black hole binaries can really coalesce 
efficiently is however a matter of debate. In a stellar en- 
vironment, it has been argued that the binary hardens only 
very slowly jBegelman et alj [l980l : iMilosavhevic & MerrittI 
l2003h . while in gaseous environments binaries may coa- 
lesce ra pidly owing to strong dynamical friction with the gas 
jMakin o & Funato 20Q4; E scala et al. 2004). In our galaxy- 
sized simulations, and even more so in the cosmological 
boxes, it is not possible to treat in detail the problem of bi- 
nary hardening, nor to directly calculate the ejection of black 
holes by gravitational recoil, or by three-body sling-shot ejec- 
tion of black holes in triple systems. Because galaxies have 
typically large central concentrations of gas we instead as- 
sume that two black hole particles merge quickly if they come 
within the spatial resolution of the simulation and their rela- 
tive speed lies below the gas sound speed. In practice, this 
means that two sink particles that fulfill these conditions are 
merged into a single BH particle, with their masses combined. 

2.3. Simulation runs 

Cosmological simulations which include quasar formation 
must model sufficiently large volumes to sample a representa- 
tive part of the universe, but also have high enough resolution 
to model the full hydrodynamics. This is a substantial chal- 
lenge, given that the brightest quasars at z 6 have a low 
space density and are believed to reside in fairly massive dark 
matter halos of mass M ^ 10" — 10'^ Mq, or even larger. At 
redshifts z = 2 — 3, quasars have a much larger space density, 
comparable to L* galaxies at z = 0. 

The strategy we choose here is model the universe with 
a periodic box of moderate size 33.75 /i^'Mpc that is ho- 
mogeneously sampled with particles, and which we simulate 
with 2 X 486-^ particles, one of the highest resolutions so far 
achieved in a full cosmological hydrodynamical calculation of 
galaxy formation. In this paper, we refer to this largest sim- 
ulation among our simulation set as the BHCosmo run. We 
will also compare it with two additional simulations which 
differ in mass and spatial resolution, and/or box size, to test 
for resolution effects. 

The fundamental numerical parameters of our simulation 
runs are listed in Table [T] where Np is the number of dark 
matter and gas particles, wdm and mgas are their masses (initial 
mass in the case of the gas). Finally, e gives the comoving 
gravitational softening length, and Zend the final redshift of the 
simulation. 

For the physical problem at hand we prefer relatively high 
resolution in order to capture the physics in high density re- 
gions appropriately, but we also need a large volume to study 
the growth of deep gravitational potentials. Our choice of 
33.75/i^'Mpc represents a compromise in this respect. While 
this box-size is too small to be evolved to redshifts lower 
than z — 1 (otherwise the fundamental mode would become 
non-linear), it is sufficiently large to provide a representative 
model for Lv,-objects at higher redshift, even though very rare 
systems in the exponential tail of the mass function will not 
be sampled well. We also note that the choice of box-size and 
particle number in the BHCosmo run is such that the physi- 
cal resolution at z ^ 6 is comparable to that in some of our 
previous works on galaxy mergers, namely runs which used 
only ~ 10000 particles for each galaxy. In this prior work. 
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Fig. 2. — Three level zooms into the simulation region marked by the white rectangle in the z = 3 panel of Fig[T] The three panels show the gas surface density, 
color-coded by temperature. The panels show slices of thickness 5000/!^' kpc and of decreasing width (from A to C) as we zoom into the region around a black 
hole of mass, Mbh ~ 7 x lO'. The yellow circles in the bottom right panel show the black holes in the region, with a symbol size that is related to the BH mass 
as in Fig.[T] 



we have shown that despite the low resolution the results for 
the black hole mass growth agreed well with those obtained 
in runs with 128 times higher resolution, and can therefore 
be considered converged with respect to this quantity. This 
overlap in resolution between our cosmological runs and our 
previous work on isolated galaxy mergers gives us confidence 
that the results of our cosmological runs are not dominated 
by resolution effects, although it is clear that this needs to be 
tested separately. We remark that as part of our previous work 
on mergers we have run a suite of several hundred galaxy 
merger simulations (Rob ertson et al.ll200 6b'). varying all the 
parameters describing star formation and feedback from su- 
pernovae and black hole growth and accretion, besides carry- 
ing out numerical resolution studies. The galaxy merger sim- 



ulations are clearly much better suited for investigating the 
full parameter space of our model, while for the cosmological 
runs we have to restrict ourselves to our default model owing 
to their much larger computational cost. 

In addition to the above considerations, the choice of box- 
size in our new simulations is also motivated by the set of 
simulations presented in Springel & Hernquist (2003b). In 
fact, the BHCosmo run would be called 'D6' in their naming 
scheme. Being able to directly refer to their runs simplifies 
the comparison of the physical properties of simulations with 
and without black holes, e.g. with respect to the star formation 
history. 

3. RESULTS 
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Fig. 3. — Stellar density in a three level zoom into the same regions shown in Figure|2] and which is marked by a white rectangle in the z = 3 panel of Fig[T] 
In the most zoomed-in panel C, a stellar disk and a small bulge component can be seen for the central object. 



3.1. Visualization of the structure and black hole growth 

In Figure [T] we show slices through the BHCosmo simula- 
tion at a range of redshifts in order to visualize the evolution 
of the baryonic density field and the growth of black holes. 
The slice has a thickness of 5/i^'Mpc and shows the full box 
of size 33.75/z~'Mpc on a side. In each panel, the projected 
gas density field is color-coded according to the gas temper- 
ature, with the brightness of each pixel being proportional to 
the logarithm of the gas surface density. Circles of different 
size are drawn to mark the locations of BHs of different mass, 
as labelled. 

The images show that black holes emerge in halos starting 
at high redshift (as early as z ~ 12) and subsequently grow by 
gas accretion, driven by gas inflows that accompany the hi- 
erarchical build-up of ever larger halos through merging. As 
the simulation evolves, the number of black holes rapidly in- 



creases and larger halos host increasingly more massive black 
holes. We note that the particular slice of the box shown in 
Figure [T] does not contain the largest black hole in the simu- 
lation's volume, which turns out to be located in the highest 
density region in the simulation. We will discuss this region 
separately in later sections. 

By plotting the density field color-coded by temperature we 
can also see traces of heating effects from strong gas outflows, 
which are caused by black hole feedback and galactic winds 
from star formation. In our numerical mode, quasar-driven 
outflows occur once a central BH gas grown so much that its 
energy feedback in accretion phases is able to transfer suffi- 
cient energy to the remaining gas to unbind and expel part of 
it from the galaxy potential as a wind. This then terminates 
further strong growth of the BH, which is key factor in estab- 
lishing a self-regulated nature of the growth of black holes in 
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our model, as we discussed in detail in our previous work on 
galaxy merger simulations. 

To illustrate how well the mass resolution of the simulation 
captures details of galaxy formation sites, we zoom in onto a 
region of 6/;^'Mpc on a side at z = 3 (as indicated by the box 
drawn in Fig[TJ, showing the gas density in Figure ID and the 
stellar density of the same region in Figure|3] The middle pan- 
els (labeled B) in both figures show a second zoom-level into 
a region of 2/z~'Mpc on a side (indicated by the square in the 
A panels). Finally, the top panels show a further enlargement 
by a factor of 8. 

In accordance with our seeding procedure, black holes are 
located in the highest density regions. The more massive 
ones are found within the largest halos, which have also un- 
dergone more prominent star formation. Such a correspon- 
dence is expected in our model since large-scale gas inflows 
into the centers of halos lead both to star formation (and star- 
bursts) as well as to nuclear black hole growth. In the highest 
level zoom of Figures|2]and|3] the central galaxy, with an ex- 
tent of ^ 50/i^'kpc (comoving), has a very rough disk-like 
morphology with a central stellar bulge. Nevertheless, it is 
clear that our cosmological simulations in general still have 
too low resolution for properly resolving galaxy morpholo- 
gies. We also note that producing disk galaxies with the right 
size and abundance in cosmological hydrodynamical simula- 
tions is an essentially unsolved problem, and the outlook for 
obtaining a solution to this long-standing challenge has only 
slightly improved by recent works on disk galaxy formation 
jRobertson et al.ll2004l; iGovernato et al.ll2007l; lOkamoto et all 
l2007h . 

3.2. The evolution of the global black hole mass density 

The black hole mass density at the present epoch is esti- 
mated from direct measurements of black hole masses in local 
galaxies (to establish, e.g., the Mbh — relationship), com- 
bined with a suitable in tegration over the g alaxy luminosity 
function (Fabian 1999; Yu & Trem aine 2002; M arconi et akl 
12004) . The local value thus obtained matches that of the to- 
tal relic BH mass density estimated from a time integration 
of the luminosity output of active galactic nuclei and quasars 
at X-ray and optical wavel engths ([Soltan 1982; M arconi et al.l 
HoollShankar et al. 2004), or in a less model-dependent man- 
ner from an empirical determination of the bolometric quasar 
luminosity function (Hopkins et al.ll20"07bh . 

Figure|4]shows our simulation prediction for the global den- 
sity pbh and its evolution with redshift {thick black line). We 
find that the normalization of the black hole mass density 
is in agreement with the observational estimate of Pbh(O) — 
4.6;!:4 X IO^/i^ vMq Mpc"^^ of iMarconi et"!!] (l2004h and its 
extrapolation to z 3, derived by exploiting hard X-ray and 
optically selected AGNs and quasars. The grey area in Fig- 
ure |4] shows the region delimited by obse rvational constraints 
examined in detail i n the literature (e.g ISalucci et al.l 119991: 
IMarconi et all2004l; IShankar et al.ll2004l) . The agreement we 
see here is very reassuring. It shows that our black hole accre- 
tion and feedback model in the BHCosmo run is adequate and 
provides a realistic account for the dominant mode of global 
black hole growth in our universe. 

The black hole mass density evolves rapidly at high red- 
shift, increasing by four orders of magnitude between z ~ 10 
to z ~ 3. While below this redshift, there is some further 
growth, it only accounts for roughly a doubling of the BH den- 
sity down to z = 1 . This is corroborated by the bottom panel 
of Figure m where we plot the history of the global black hole 
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Fig. 4. — Top panel: The global black hole mass density evolution in 
the BHCosmo/D5 run, shown by a solid line. The the star symbols and the 
dashed line give the corresponding stellar mass density evolution, multiplied 
by 0.0007 for easier comparison. The grey dotted line shows the stellar 
mass density in the D5 simulation which did not include black hole accre- 
tion and associated quasar feedback (Sprinsel & Hernauist 2003a). The thin 
dashed line shows the results from the lower resolution box, D4 described 
in Table [T] Different colors simply indicate the different redshifts consis- 
tent with the scheme used in other figures. The shade d grey triangle indi- 
cates observational constraints taken from the literature dMarconi et al.l2004t 
IShankar et al. 2004). Bottom Panel: The global history of the black hole ac- 
cretion rate (solid line) and star forination rate (dot-dashed line with stellar 
symbols) densities. The SFR is rescaled by 10^^ for graphical clarity. In 
addition, we show the SFR history in the D5 simulations without black holes 
(grey dotted line). Most of the black hole and stellar mass is assembled by 
z 2 — 3, but the peak in the BHAR density function is far more pronounced 
than that of the SFR density. 



accretion rate (BHAR) density in our BHCosmo run. It is 
evident that the black hole mass growth occurs mainly by ac- 
cretion above z ~ 3. Indeed, the BHAR density rises steadily 
at early times to a peak at z ~ 3, where it reaches a value 
^ lO^'^M© yr^' Mpc^"'. Below this redshift, it drops rather 
rapidly, becoming more than an order of magnitude lower by 
z = 1. We note that the 'spike' seen in the BHAR density at 
z ~ 6.5 is caused by a single object in the box, namely the 
rapid formation of the most massive black hole in the simu- 
lation at this epoch. We will discuss the history of the corre- 
sponding halo and its embedded black hole in more detail in 
§5. 

3.3. Black hole mass function and accretion rate function 

In Figure|5] we plot the black hole mass function at a num- 
ber of different redshifts. We find that the final black hole 
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Fig. 5. — The evolution of the black hole mass function in the BHCosmo 
run. The different colors indicate different redshifts as labeled in the figure. 
For comparison, the dark grey shaded region shows the black hole mass func- 
tion derived from local galaxies (Marconi et al. 2004). The light gray area 
adds th e constraint from the integration of the hard X-ray luminosity func- 
tion (S hankaret"ai] [2004: Merlonil l2004l : lMarconi et al.i .2004). The largest 
uncertainties are at the high and low mass end. The simulation results are in 
good agreement with the observed mass function at the high mass end, and 
in reasonable agreement at intermediate masses. 
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log (m/rhEdJ 

Fig. 6. — The time evolution of the accretion rate distribution as a function 
of the Eddington ratio, for the BHCosmo run. The different colors denote our 
measurements at different redshifts, as indicated in the legend. For the z = 1 
distribution function, we separately show three components corresponding to 
different regions of the black hole mass function. In particular, the dotted, 
dot-dashed, and dot-dot-dashed lines give the separate contributions from the 
three different mass bins Mbh > 10*^ 
Mbh/Mq < 10^, respectively. 



Mq, lO' <Mbh/Mq < 10** and 10^ < 



mass function in our simulation (for z = 1) is is quite good 
agreement with the one measured locally, especially on the 
high-mass side. The z = constraint is indic ated by the 
dark grey area taken from the compilation of M arconi et alj 
(l2004i) . which is based on a combin a tion of different obser- 
vational data dKochanek et all 120011: iNakamura eTaTI 120031: 



iBernardi et all 120031: ISheth et all l2003h . There is a small 
deficit in our model at intermediate BH masses, but note that 
this will be filled in at least partly by the expected residual 
growth from z = 1 to z = 0. The light grey area adds an ad- 
ditional constraint for the contribution of relic AGN, derived 
from an integration o f the the hard X-ray luminosity function 
dShankar et alj2004l) . (Note that in this latter case the normal- 
ization of the mass function depends on the value assumed for 
the radiative efficiency when converting from the luminosity 
function to the mass function.) 

As expected in a hierarchical formation scenario, the high 
mass end of the mass function shifts to larger masses with 
redshift. However, it is interesting that the mass function for 
high masses grows rather rapidly at early times relative to the 
low mass end. On the other hand, for redshifts z < 2, the 
high mass end is virtually fully assembled while the BH mass 
function continues to grow for low to intermediate masses, 
leading effectively to a steepening of the mass function for 
z ^ 3. This appears consistent with the emerging observa- 
tional picture of the evolution of the supermassive black hole 
population according to which massive BHs (Mbh ^ 10"^ Mq) 
are assembled early and are then likely undergoing compar- 
ativel y passive evo l ution in the centers of large spheroids 
(_ShieldsetalJl2003l: IVestergaardll2004t lAdelberger & Steidell 
I2OO5I) . More generally, this phenomenon is described by the 
idea of an 'anti-hierarchical' black hole growth, or equiva- 
lently a 'downsizin g' of black hole activity (Merloni 20041 
iMarconi et al.ll2004h . which is derived from constraints on the 
accretion history from X-ray luminosity functions. 

Finally, in Figure |6l we show the evolution of the cor- 
responding accretion rate distribution function for the black 
holes in our simulation, expressed in units of the Eddington 
rate. This function is strongly peaked at a few percent below 
the critical Eddington value, with most black holes accreting 
at 10^^ < m/»jEdd ^ 1 at redshifts z^6. The distribution 
becomes wider and develops a small amplitude tail at low Ed- 
dington rates for 6 < z < 3. Below z ^ 3 the peak of the ac- 
cretion rate function shifts to m/rh-EM ^ 10^^ — lO^"* and an 
increasingly large population of sources is present accreting 
down to 10^^ Eddington. At z 1 the distribution function is 
extremely wide and dominated by sub-Eddington (to severely 
sub-Eddington) sources, in terms of number. While hence the 
number of sources accreting at sub-Eddington rates increases 
sharply with decreasing redshift, the quasar-like population, 
i.e. the sources accreting close to Eddington, decreases with 
decreasing redshift, with a maximum at z 3. For z ~ 1, we 
also include in Figure |6] measurements of the separate con- 
tributions to the accretion rate function from three different 
black hole mass ranges. The dotted line gives the Eddington 
ratio distribution for Mbh > 10*^ Mq, and the dash-dotted and 
dot-dot-dashed lines are for masses 10^ < Mbh/Mq < 10** 
and 10'' < Mbh/Mq < 10^, respectively. 

Taken together, the evolution of the black hole mass and 
accretion rate functions implies that most massive black holes 
assemble early and do so at close to their critical rate. At low 
redshift, progressively lower luminosities and lower mass sys- 
tems start dominating the black hole activity. This is in accord 
with recent results from studies of the hard X-ray luminosity 
function of quasars and AGN, e.g. typical Seyfert-like objects 
with Mbh 10' — 10^ Mq acc reting at a few percent of the 
Eddington rate (e.g. lUeda et al.i,2003; Hasinger et al. .2005) . 
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Fig. 7. — An illustration of a large group in the BHCosmo run at z = 1. Left panel: the gas density distribution. Right panel: the stellar distribution. We have 
run a sub-group finder to identify all the systems (galaxies) within each halo and analyze their properties accordingly. The stellar system shown in yellow on the 
right hand side is the main galaxy within this large group, while the satellite galaxies are shown in grey. The images are 400 kpc on a side. 



3.4. Comparison of the histories of black hole growth and 
star formation 

In Figure m we have already shown and discussed the evo- 
lution of the black hole mass and accretion rate densities. In 
the same Figure, we show corresponding results for the evo- 
lution of the stellar mass density, p*, and the star formation 
rate (SFR) density, with their normalizations rescaled by fac- 
tors of 7 X 10^^ and 10^^, respectively, for plotting purposes. 
F or comparison, we also include results for the D5 simulation 
of lSpringel & HernquistI (l2003bh which did not include black 
holes and any associated accretion or feedback processes (dot- 
ted grey lines). 

We can see that far exceeds pbh at all redshifts, with 
evolving less strongly with redshift than p^n for z 3. At 
early times, pbh rises more rapidly than the star formation 
density, while it tracks its evolution below this redshift. If 
we parameterize the ratio of Pbh/p* by an evolutionary factor 
(1 + z)" we find that its evolution is approximately given by 

r 03(1 +z)-^ ifz>3.5 
Pbh/p* ^ I , (1) 

I '?!>i[(l+z)]"°-^ ifz<3.5 

with (/)3 = 5 X 10~^ and 4>\ — 2.2 x lO^"*. Accordingly, up 
to z ~ 3, the evolution of the star formation rate density is 
considerably shallower than that of the black hole accretion 
rate density. Below this redshift, the BHAR and SFR densities 
closely track each other As a result, the BHAR density has 
much more pronounced peak, which we find to lie at slightly 
lower redshift than the peak of the SFR. Parameterizing the 
evolution of the BHAR and SFR density ratio with a power 
law in (1 +z) leads to 

r P7.{\+z)-'^ ifz>3 
Pbh/p* ~ < , (2) 

I ifz<3 



with 03 0.2 and 0i 10^1 

Note that Equations ([U and (|2|l are only meant to provide 
approximate scalings for our results from the simulations. 
The important point we want to emphasize is that our results 
imply a different and much stronger evolution of the black 
hole mass and accretion rate densities at high redshift rela- 
tive to the stellar density and star formation rate density. The 
black hole mass density tends to be assembled later than the 
stellar mass, despite the growth of (a small number of) very 
massive BHs already at high redshift. However, for z 3 and 
below, our models predict that the black hole mass and stel- 
lar mass densities, as well as the BHAR and SFR densities, 
closely track each other allowing only for small amounts of 
relative evolution between the two. 

Another important point from Figure |4]is that the feedback 
we associate with black hole accretion does not significantly 
affect the global assembly of stellar mass. The peak of the 
SFR density is unaffected by the inclusion of BH feedback, 
but the drop in SFR density (z < 3) is slightly more abrupt in 
the simulations with black holes. This effect becomes more 
pronounced at z ~ 1, the final redshift for our simulation. At 
still lower redshift, we expect that BH feedback will become 
important in regulating the cooling and star formation in very 
massive halos. This is then ascribed not to quasar growth but 
rather to a 'radio mode' . We explore this different feedback 
channel in a companion paper by Siiacki et al. (2007). 

4. THE Mbh - cr AND Mbh - M, RELATIONS 

4.1. Identification of groups and subgroups 

As a prerequisite for being able to study correlations be- 
tween black hole and host galaxy properties in our simulation 
we first need to apply a suitable group finding algorithm that 
reliably identifies the stellar mass associated with the different 
galaxies. Note that especially the more massive halos identi- 
fied by our basic friends-of-friends grouping algorithm used 
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Fig. 8. — The evolution of the Mgn ^ f relation in the BHCosmo simulation. The masses of BHs and the projected stellar velocity dispersions within the half 
mass radius (Re) have been measured in our simulated galaxies and are plotted at z = 1, 2, 3. 4, 5, and 6.5. We compared our measurements at all redshifts with 
the best-fit to the local Mbh — f relation of Tremaine et al. 1 2002, hereafter T02), which is shown as a thick gray line. Linear regression fits to our simulated BHs 
are shown by solid fines at each redshift, with \ — a errors indicated by the hatched regions. For ease of comparison, the dotted lines in each panel show the best 
fit relations at all redshifts. The points are color-coded according to their accretion rates in units of Eddington, as indicated in the color bar at the top right hand 
comer of the figure. 



for finding virialized objects often contain a number of galax- 
ies. This is illustrated in Figure |71 where we show a large 
cluster-sized group selected from the z = 1 output of the BH- 
Cosmo simulation. The panel on the left shows the gas density 
distribution in this large group, while the panel on the right 
hand side displays the stellar distribution. It is evident that 
the group contains several, gravitationally bound galaxies. 

Our sub-group finder identifies all the galaxies within each 
group. Our method to identify galaxies within a given group is 
based on a variant of the SUB FIND algorithm (Spring el et al.l 
I2OOII) . We first compute an adaptively smoothed baryonic 
density field for all stars and gas particles, allowing us to 
robustly identify centers of individual galaxies. We find the 
extent of these galaxies by processing the gas and star parti- 
cles in the order of declining density, adding particles one by 
one to the galaxies attached to the galaxy to which its near- 



est de nser neighbor already belongs (see also'Nagamin e et aP 
l2004b . Note we are interested in the stellar and gas content 
of galaxies and associating the gaseous component to galax- 
ies, as they typically they contain very dense star-forming gas, 
makes the method very robust in finding galaxies. This allows 
us to compute physical properties such as stellar mass, star 
formation rate, stellar velocity dispersion, metallicity, black 
hole mass and BH accretion rate separately for each galaxy. 
As an illustration, Figure|7]shows the stellar distribution asso- 
ciated with the largest galaxy in the group (which also hosts 
the most massive BH in the group) in yellow. 

For each galaxy/subgroup that contains stars and a black 
hole, we calculate the projected (spherically averaged) half- 
mass effective radius, R^, and the mass-weighted stellar ve- 
locity dispersion a within R^. Note however that in order to 
make the measurements of a and Re somewhat accurate we 
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Fig. 9. — The evolution of the Mbh — 1^* relation in the BHCosmo simulation. The masses of BHs and the con'esponding stellar mass have been measured in 
our sim ulated galaxies and are plotted a.tz= 1 , 2, 3, 4, 5, and 6.5. They are compared with the best-fit power law to the local Mbh — M, relation by Haring & Ri^ 
I20Q4I) . shown by the thick gray line. Our fits to the simulation results at each redshift are shown as solid black lines. The relation is tight with small scatter 
(1 — cr uncertainty ranges are plotted as hatched regions but are hardly visible at low redshift). As in Figure[8] the points are color coded by accretion rate. The 
dotted lines show in each panel the best fit relation at the other redshifts. 



only consider those objects that contain more than 100 stars 
particles within Rg. This determination of a closely resembles 
the procedure fo r measuring the veloci ty dispersion from ob- 
servational data jGebhardt et alj|2000l) . allowing for a direct 
comparison. 

4.2. The predicted Mbh — cr ond Mbh — M* relationships and 
their evolution 

Figures [8] and |9] plot the Mbh — f and Mbh — M* relations, 
respectively, for our simulated galaxies at redshifts z— 1,2, 
3, 4, 5 and 6.5 (from top left to bottom right). Each mea- 
surement is color-coded according to the accretion rate of the 
corresponding black hole. We find a strong power-law cor- 
relation between both the velocity dispersion and the stel- 
lar mass with the black hole mass. Furthermore, at z = 1, 
these correlations agree very well with the ones observed at 
the present epoch dFerrarese & Merritj 120001: iGebhardt et alj 



l2000tlTremaine et al .ll2002l: iHiiring & Rixll2004l) over a large 
dynamic range. This is a remarkable confirmation of the basic 
merger-driven scenario for self -regulated BH growth that we 
previously explored in isolated high-resolution galaxy merger 
simulations. Based on the same model, and without any fine- 
tuning or change of the model parameters, we here obtain a 
good match to the observed Mbh — c relation directly from 
simulations that start from cosmological initial conditions and 
self-consistently account for the hierarchical build up galaxies 
in a ACDM universe. 

4.2.1. The evolution of Mbh — (7 

At higher redshifts, and in particular for z > 3, our results 
for the Mbh ^ cr relation shows some degree of evolution rel- 
ative to the local relations. To compare with the local data as 
a function of time, we fit the Mbh ~ c relation at each redshift 
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TABLE 2 

Parameters of best-fit Mbh - ^ relations 



z slope (a) normalization, (b) scatter A slope (aj) 



1 


3.95±0.10 


8.29 ±0.03 


0.10 


3.9 


2 


4.01 ±0.15 


8.42 ±0.04 


0.16 


4.1 


3 


4.21 ±0.22 


8.32 ±0.07 


0.23 


5.2 


4 


4.54 ±0.35 


8.17±0.10 


0.36 


6.7 


5 


3.37 ±0.45 


7.70±0.13 


0.47 




6=; 


.. 3.26±0.85 


7.50 ±0.25 


0.89 





with a simple power-law of the form 

log (^]=a log f ^) +b. (3) 

^\MqJ ^V200kms-i / ^ ^ 

The best-fit relations thus obtained are shown in Figure [8] as 
solid black lines, with the dispersion indicated by hatched re- 
gions. The dotted lines are the best-fit relations for all red- 
shifts combined. We compare with the observed relation as 
determined by T02, which is described by a slope a — 4.02, a 
normalization b — %.2 and a dispersion A ~ 0.25 — 0.30 (grey 
thick line in Fig.[8]l. 

The constants a and b for our best-fit relations and their dis- 
persions are tabulated in Table|2]for different redshifts. Note 
that we here do not attempt to assess the statistical signifi- 
cance of the correlations in detail as the sources of system- 
atic errors in the numerical measurements of a and Mbh can- 
not be easily quantified in the simulations. In particular, the 
cosmological simulations cannot determine the morphologi- 
cal properties of galaxies and therefore do not provide a direct 
measure of spheroid masses or their velocity dispersions. Our 
fitting procedure is merely intended to provide a first char- 
acterization of the overall evolution of the slope and normal- 
ization of the relations in the simulation model. As Figure |8] 
and Table|2]indicate, the Mbh ~ f relation predicted from our 
simulation is consistent with a slope ~ 4 at low redshift, as ob- 
served. At z = 3 — 4, the slope appears to be slightly steeper 
and at z = 5 — 6 slightly shallower, but the small number of 
systems at z ~ 6 makes the latter trend uncertain. 

Inspection of Figure |8] shows a qualitative trend whereby 
the larger systems with high a ^ 150 km s^' (which are also 
those that are best resolved in our simulation) appear to pop- 
ulate the high mass end of the Mbh ~ f relation already at 
z > 2 — 3. The lower mass end of the Mbh — f relation is 
then increasingly filled in towards z 1. To illustrate this 
trend more explicitly, we perform a fit of the relation using 
only those systems with a > 150 km s^' ; this is shown by 
the dashed line in Figure |8] The slope Oj obtained for these 
'high' (T black holes is typically steeper than for the full rela- 
tion (our measured values for Oj are listed in Table |2|i. At a 
fixed and relatively high a, the mean Mbh is larger at z 3 — 4 
than at z '--^ 1. This trend suggests that black hole growth 
predates the final growth of the spheroid potential at scales 
(7 ^ 170kms^', which is consistent with the recent mea- 
surements of t he Mbh ~ o ' relation in a sample of Seyferts 
at z 0.35 by IWoo et al.l (12006) . This relatively more ef- 
ficient black hole growth at a fixed c at high redshift may 
thus be responsible for a mild change in the high mass end 
of the Mbh — cr relation. In the remainder of this section we 
will discuss possible additional dependences in the Mbh — f 
jHopkins et al. 2007a) relation which may produces a small 
change in the slope of the relation with time. 




M . (Mo) 



Fig. 10. — The stellar velocity dispersion a versus stellar M* at z = 1,2, 
3, 4, 5, 6.5 (i ndic ated by different colors from blue to pink, the same ones as 
used in Figs. 14161 . The best-fit power-law to the trend is shown with a dotted 
line at each redshift, and with a solid line at z= 1 . The dispersion a at fixed 
M, increases with increasing redshift, which can be interpreted as a weak 
evolution in the Faber & Jackson relation. 



Note also that at z > 4 black holes are more likely to accrete 
close to their critical Eddington rates, as we showed earlier. 
In our models, the Mbh — cr relation is a natu ral consequence 
of the self -regulated growth of black holes (iPi Matteo et alj 
l2005h . where a black hole grows until its released energy 
is sufficient to expel the gas in its surroundings in a quasar 
driven wind, which then terminates nuclear accretion. For this 
reason, we expect the relation to show more scatter at times 
when most systems are still actively growing and accreting 
close to their Eddington values (see Table|2l). This is expected 
as the primary path for assembly BH mass is via accretion 
during major mergers so that the relations converge and get 
increasingly tighter as galaxies undergo major mergers and 
continue to merge. 

The overall normalization of the Mbh — o" relation evolves 
weakly with redshift when compared to the T02 result for 
the present epoch. Below z 3, there is virtually no evo- 
lution (bearing in mind the weak trend discussed above for 
(7 > 150 kms^' ), while we find a weak drift of the normal- 
ization beyond z 3. A rough parameterization of this evolu- 
tion in the normalization of the Mbh — c relation is given by 
Mbh/ct"^ oc (1 + z)" with a = -0.2. 

These results for an overall modest evolution in the normal- 
ization are in qualitative agreement with tho se from an anal- 
ysis of i solated galaxy merger simulations bv lRobertson et al] 
(l2006b!) . We note that in the latter study a redshift scaling 
of galaxy properties and a specific set of structural properties 
of the merging galaxies had to be assumed, which introduced 
an important systematic uncertainty. While having lower nu- 
merical resolution per merger, the simulation we analyze here 
eliminates this caveat by providing a fully self-consistent cos- 
mological history for the formation and evolution of galax- 
ies and their black holes (albeit at a lower spatial resolu- 
tio n). This provides an imp ortant confirmation of the analysis 
of [Robertson et al] (l2006bl) . who however did not find evi- 
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Fig. 1 1 . — The redshift evolution of the projected half mass radius Re (top 
panel), cold gas fraction /ga.s (middle panel) and characteristic black hole to 
stellar mass ratio {M^n/Mf) (bottom panel), for systems with M^, > 3 X 
lO'^M© (dashed lines) orM* > 6 X IO'^Mq (solid) in the BHCosmo run. 
These threshol d values were chosen to compare with the observed evolution 
determined bv lTruiillo et al.l I200flh . The increase in cold gas content in high 
redshift progenitor hosts, and the trend to more compact systems with an 
increasing ratio of M^w/M^ at a fix ed s tellar mass is con sistent with the 
recent results of Trujillo et al. ( 20061) and lPeng et alj j2006l) . as well as the 
BHFP iHopkins et al..2007ai) . 



dence for an evolution of the slope at the high mass end of the 
A^BH — 0" relation. 

4.2.2. The evolution of Mbh - 

In Figure |9] we show the Mbh — relation from the sim- 
ulations at z = 1, 2, 3, 4, 5 and 6. 5, alongsi de the local ob- 
servational relation determined bv iHaring & Rix, (2004 , thick 
grey lines). Our best-fit relation at each redshift is shown by a 
solid line while the dotted lines in each panel show results for 
the other redshifts. Table[3]gives the slope c, normalization d, 
and dispersion A,„ for our best-fit relations of the form 
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1011 M, 
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As before, our fitted values for c and d are intended to indicate 
general trends in the evolution rather than to be used as sta- 
tistically rigorous measurements. The observed relationship 
(|Haring & Rix 2004) has a slope c = 1.12 and normalization 
£/ = 8.2. 

Overall there appears to be only limited evolution in the 
Mbh — M, relation, but there is a slight steepening at z = 2 — 4. 
To highlight this trend, we restrict our fits to the high mass end 
withM* >5 X 1O1"M0 (dashed line inFig.|9ll. In this range, 
the relation is significantly steeper, implying slopes c, ~ 1 .9 
at z = 3 — 4 and ~ 1 .5 at z ~ 2. This is more significant than 
the evolution found in the slope of Mbh — c, and implies that 
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there is some evolution in the ratio of black hole mass to stel- 
lar mass relative to the local observations. More precisely, 
systems with M» ^ lO"' Mq have larger black hole masses at 
fixed M* than at z = 1, where the ratio is in good agreement 
with the relation observed at the present epoch. This trend of 
an increasing Mbh/M* ratio as a function of redshift appears 
consistent with the recent measurements of high redsh ift (up 
to_z_~ 3.5) BH masses and h ost luminositi es by Pen g et alj 
( I2OO6I) . as well as the BHFP (Hopk ins et ani2007ai) . We will 
further analyze this effect in § 14.31 

4.3. Evolution of a, gas fraction, and Re, at fixed stellar host 

mass 

We now analyze some of the physical properties of the host 
galaxies and their evolution with redshift to investigate the 
physical origin for the trends we have found in the Mbh ^ cr 
and Mbh — M^, relations. Figure [TOl shows the stellar veloc- 
ity dispersion a versus the stellar mass M* for each galaxy 
as a function of redshift. The dotted lines and the solid line 
(for z = 1) show our best-fit relations. At a fixed M*, the ve- 
locity dispersion a increases with increasing redshift. This 
is consistent with the r esults from the merger remnants in 
iRobertson et al.l (l2006bl) and reflects changes in the structural 
properties of stellar spheroids, which become smaller towards 
higher redshift. 

In Figure [m we show the projected half mass radius R^, 
the cold gas fraction /ga,s (a proxy for the disk gas fraction), 
and the (Mbh/M,^) ratio as a function of redshift. We show 
results for systems with M* > 3 x lO"' Mq (dashed lines) and 
M* > 6 X IQi" Mq (soUd fines). We find higher cold gas frac- 
tions in higher redshift systems with fgas (x (1 +z)''^. The 
increased gas content, dissipation rate in high redshift progen- 
itor is consistent with them being more centrally concentrated, 
as shown by the trend of decreasing R^ and increasing cr with 
increasing redshift, leading to larger ratios Mbh/M* at high 
redshifts, with Mbh/M* ~ (1 +z)°'^, at least over these ranges 
of M* . Larger values of a at fixed M* at high redshift are also 
consistent with an inverse overall trend (e.g.; (1 +z)^°^) in 
the evolution of the Mbh — f relation. 

The above results are consistent with the general picture 
for a merger-driven quasar growth we have developed based 
on our galaxy merger simulations (Di Matteo et al.' '2005t 
Springel et al. 2005b; Robertson et al. 2006b; Hopkins et al] 
2005. i2006al I2007al) , which also implies that the Mbh - cr 
and Mbh —M* relations are connected to each other. This con- 
nection is most clearly demonstrated by the ' black hole fun- 
damental plane' (BHFP) relation discussed bv lHopkins et alj 
(2007a), which arises from the joint formation process of 
spheroids and massive black holes in mergers. Our cosmolog- 
ical simulations also should agree with the BHFP if most of 
the BH growth is associated with mergers. In Figure [121 we 
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Fig. 12. — Simulation prediction for two projections of the 'B lack Hole Fundamenta l Plane' (BHFP) relation at z = 1, in terms of a and M, (left panel), and a 
and Re (right panel). We compare with the best-tit relations from lHopkins et al] j2007at) . shown as dotted lines. The simulation agrees well with the conjecture of 
a BHFP, which confirms the overall trends we have found in the Mbh — f and ^BH — ^* relations. This likely owes to the scatter in the measurements of both a 
and in particular Re which are noisy measurements in the simulations. The colors indicate accretion rate values (as in previous figures and indicated in the color 
bar). 



show our simulation measure ments for the BHFP a nd com- 
pare it to the the best-fit from iHopkins et aP (l2007ah . shown 
as a solid line. We find very good agreement w ith the predic- 
tions for the BHFP by iHopkins et al] (l2007al) . in both of its 
formulations, i.e. for Re and a or and a, although the for- 
mer shows a larger scatter than the latter, owing to the noisy 
measurements of in the cosmological simulation. 

5. THE FIRST AND THE MOST MASSIVE BLACK HOLES 

In the previous sections we have discussed the predictions 
from our simulation for the global history of black hole mass 
assembly in galaxies from the high redshift Universe to to- 
day. We have compared the predictions for the evolution of 
the black hole mass and accretion rate density to the history 
of the star formation rate density and discussed the growth 
of the black hole mass function. We have also discussed the 
A^BH — f and Mbh — * relationships as a function of time, 
and examined which physical properties drive their cosmo- 
logical evolution. 

However, our simulation methodology not only allows us 
to make statistical statements about the black hole population. 
Rather, it can also be used to study the detailed growth history 
of individual black holes, from the moment they are seeded to 
today, which provides a particularly powerful way to follow 
the evolution of black holes over cosmic time. The gas which 
fuels black hole activity ultimately has its origins in the inter- 
galactic medium, draining along filaments into forming galax- 
ies. Because of this, BH radiative histories are directly linked 
to the formation of large-scale structure in the Universe, from 
supercluster scales down to the immediate environment of 
host galaxies. Being able to follow these large-scale processes 
and their impact on individual black holes self-consistently is 
the key to a qualitatively better level of understanding. For 



example, we have the tools to examine what turns a small 
black hole into a supermassive one at z = 6, or whether some 
BHs grow hardly at all after this initial phase (as we shall see 
below in some examples). We can ask how quasar lifetimes 
are related to their clustering, how important BH mergers ver- 
sus gas accretion are in the cosmological growth of BH mass, 
what BH light curves look like in detail over a Hubble time, 
and how specific outbursts correlate with galaxy and cluster 
merger and accretion events. 

In future work we will return to address suc h ques tions 
in more detail; for example in IColberg et al.l ("2007) we 
present an extensive study of BH environments as the universe 
evolves. For now we will show that there are some striking in- 
ferences to be drawn from following the detailed histories of 
even a few black holes. As an example we choose to focus on 
the formation and fate of the first large black holes that form 
in our simulation, as well as the properties of their hosts and 
their descendants. This allows us to identify the prerequisite 
conditions for the growth of supermassive black holes already 
in the early universe. 

Observations of luminous SDSS quasars at redshifts as high 
as z ^ 6 (Fan et al. 2003) present a number of challenges for 
models of high redshift quasar and galaxy formation. Their 
low space density suggests that they reside in the rarest dark 
matter density peaks at this early epoch, yet the apparent 
lack of companion galaxies in the field has been used to 
argue that these quasars reside in far more common halos 
dCarilh et al . 2004; Willott et al. 2005). Extensive follow-up 
observations of the highest redshift (z — 6.42) quasar, SDSS 
J11148-H525 1, are starting to co nstrain the properties of its 
host galaxy (IWalter et al.l 12004 ). CO observations indicate 
a relatively small stellar spheroid, however, the existence of 
CO, iron and carbon emission indicates a heavily enriched 
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Fig. 13. — Individual mass assembly and accretion rate histories for the six most massive black holes and two intermediate mass BHs (chosen randomly) in 
the BHCosmo run. The top panel shows the black hole mass as a function of redshift while the bottom four panels give the detailed accretion rate history for the 
black holes with the corresponding colors in the top panel. The first supermassive black hole forming in the BHCosmo run, and the most massive one at the end 
of the run, are shown with thicker lines in pink and blue, respectively. From the bottom panels we see that phases of high Eddington accretion occur at different 
times of the history of different black holes. 



ISM and vigorous star formation (iMaiolino et al.ll2005l) . In 
any case, rapid growth is clearly required to produce large su- 
permassive black holes and spheroids in the ~ 800 million 
years available to z ~ 6. Unsurprisingly, the type of object 
associated with the 'first quasar' is therefore still a matter of 
debate. 

The relatively small volume of our simulation prevents us 
from directly addressing the problem of the 'first quasars', 
simply because their space density is so low. Properly iden- 
tifying one of the rare hos t systems requ ires at the very least 
box-sizes of 500/i-'Mpc (Springel et alfcOOSdl) . or more ad- 
equately 1/z^'Gpc dLi et al. 2007) . Nevertheless, our simula- 
tion allows us to examine the nature of the environments and 
the hosts where exponential growth of BHs can first occur at 
early times. This has some relevance to the problem of where 
the first quasars are likely to form. 



In order to retrieve the required information for all the black 
holes in the simulation we have constructed full merger trees 
for each of them. This allows us to track the growth and ac- 
cretion history of individual black holes of different masses, 
make detailed lightcurves from their accretion histories, and 
study the properties of their host galaxies. 

5.1. Individual BH mass and accretion rate as function of z 

In Figure [13] we show the accretion histories of a number 
of black holes in the BHCosmo run (top panel) and some ex- 
amples of their corresponding accretion rate histories in the 
bottom four panels. The sample in the top panel consists of 
the six most massive black holes in the simulation and two 
randomly chosen ones with an intermediate mass at z = 1. 
The four bottom panels show the accretion rate in units of Ed- 
dington for four black holes with corresponding line colors. 





Fig. 14. — Time evolution of the environment around the host galaxy of the first supermassive massive black hole at z = 6 (bottom panels), and of the most 
massive BH at z = 1 (top panels) in the BHcosmo run. The location and masses of the supermassive BHs are marked by arrows of different size, as labeled. 
While the bottom system hosts the most massive black hole at high redshift, it does not end up hosting also the most massive black hole at the center of the largest 
galaxy at low redshift. Instead, the system shown in the top panels overtakes it in growth at intermediate redshifts, when it is formed in the highest density region 
in the simulation, which is a protocluster region. 



including the most massive black hole at z = 6, the most mas- 
sive and second most massive at z = 1, and finally one of the 
intermediate mass ones. Note that this is a small sub-sample 
of the several thousand black holes in the simulation. We fo- 
cus here on the evolution of the most massive ones, and in 
particular on their early f ormat ion and the fate of their descen- 
dants. In IColberg et al.l (l2007i) . we study a complete sample 
of individual BH histories as a function of environments for 
the entire simulation. 

The first interesting result shown in Figure [13] is that even 
in this small volume conditions exist that are conducive for 
exponential black hole growth, as required by the presence of 
quasars at z ~ 6 with the large observed masses. The evolu- 
tion of the mass of the largest black hole in our simulation at 
early times is shown by the thick pink line. The corresponding 
accretion rate history (in the same color) shows a rapid suc- 
cession of numerous phases of high accretion at the critical 
Eddington rate, between 5 z ^ 7.5. Below z ~ 5, the black 
hole is drastically more quiescent, except for a couple of spo- 
radic Eddington phases at around z^2. Interestingly, this first 
most massive black hole at z 6 is not the most massive one 
at z = 1. Instead, it is 'overtaken' in growth by another black 
hole shown with the blue thick line, which is the most mas- 
sive black hole at the end of the simulation. The Eddington 
growth phases for this objects start at z < 6 and extend all to 
way to z 3.5. By z ^ 4.5, the black hole mass in this system 
catches up with that of the first massive black hole, and then 
keeps growing exponentially to z ~ 4. 



In Figure [141 we show projected gas density maps of the 
host galaxies and the surrounding structures (at redshifts z = 
7.5, 4, and 1.6) for two of these black holes, the one that is 
most massive at z = 1, and the one that is most massive at 
z = 6. This provides good clues for the origin of the evo- 
lutionary difference between these two systems. Although 
the hosts of both black holes at z = 7.5 are halos of similar 
mass, ~ 10"^ M0, one of them lies at the intersection of three 
small filaments, leading to efficient gas cooling and a high star 
formation of ^ lOOMoyr^'. The host galaxy grows rapidly 
until z = 6 but then its gas supply dwindles, and at z ^ 4 it 
is overtaken by the other black hole. This latter one lies on 
a very massive filament, which grows even more vigorously 
at this later time than its host galaxy (which has a SFR of 
^ 1000 MQyr^' at z 2). In this way a large stellar spheroid 
is produced around the BH, which eventually will end up at 
the center of a rich cluster of galaxies. Our results therefore 
show that a massive black hole that is found in the cD of large 
galaxy cluster at late times was not necessarily the most mas- 
sive one at z = 6, as it has been often assumed in the literature 
(e.g. Springel et al. 2005d). Since the growth history of black 
holes is intertwined with the non-linear processes of structure 
formation, individual growth histories of BHs can be complex 
and need not preserve the rank order in a group of BHs that 
start out with similar masses. 

Further confirmation of this result from the BHCosmo sim- 
ulation come from our E6 run (see Table [TJ, which includes 
the same physics and model parameters as the BHCosmo sim- 
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Fig. 15. — Comparison of individual black hole mass histories in the BH- 
Cosmo run (blue lines, D6) and in the larger simulation volume of the E6 run 
(red Hnes). The growth of the first supermassive black holes at z 6 — 7 is 
more widespread in the larger volume of the E6. The 'catch-up' of larger 
black holes that form later and in higher density r egio ns (see text) can also be 
seen, and is comparable to the case shown in Fig. 1 131 



ulation but samples rare halos better, thanks to its larger box- 
size of 50/i~'Mpc on a side. In this larger volume, there are 
a handful of more examples of early exponential black hole 
growth between 8 ;S z ^ 6, as shown by the red black hole 
mass growth curves in Figure [15] Consistent with our pre- 
vious finding, there are again black holes that 'catch-up' in 
growth to the most massive black hole at some earlier time. 
We note that we have run this larger volume simulation with 
the specific purpose of checking the BHCosmo results for the 
first supermassive black holes by obtaining better statistics for 
large halos at z ^ 6. For this reason, this simulations was not 
evolved beyond z = 4. 

5.2. BH merger trees for the first and the most massive black 

hole 

Figures [16] and [TT] show two example black hole merger 
trees, one for the most massive black hole at z = 1, and one 
for the z = 1 descendant of the first supermassive black hole 
at z 6. All progenitor black holes that merge together to 
build up the final BH are included in the trees. The black hole 
masses along the tree are represented by different sizes of the 
circles, while the color encodes their accretion rate in units of 
the Eddington rate. In both figures, the top-left panel shows 
the full BH merger tree, whereas the top-right and bottom two 
panels split the tree according to accretion rate of the BHs. 

Inspection of these two BH trees immediately reveals a 
complex and rich merger history for the most massive black 
holes at z = 1 (Fig.[T6]l, as opposed to the comparatively iso- 
lated evolution of the most massive object at z = 6, which 
only features three major branches for its tree (Fig.[T7b. The 
large majority of the mass in the latter case is built up by early 
phases of critical gas accretion at z ^ 6. On the other hand, in 
Figure [16] Eddington accretion phases are spread out along 
many branches of the tree and cover a much larger range of 



redshifts. At redshifts z < 3, the progenitors are preferentially 
accreting at low Eddington rates, so that the residual growth 
of the final black hole mass increasingly occurs from "dry" 
mergers with other black holes, consistent with our earlier re- 
sults for the global black hole growth. 

6. SUMMARY AND DISCUSSION 

In this paper we have investigated the coupled formation 
and evolution of black holes and galaxies using state-of-the- 
art cosmological hydrodynamic simulations of the ACDM 
model. For the first time, we have incorporated black 
hole growth and associated feedback from quasar activ- 
ity self-consistently in cosmological hydrodynamic simula- 
tions. Our approach has been based on the methodology re- 
cently developed and tested in simulations of galaxy mergers 
(iDi Matteo et al.1 . 2005: S pringel et a l. 2005b). In this initial 
paper we have focused on investigating the model predictions 
for (i) the global history of black hole mass assembly and its 
relation to the history of star formation, from high redshift to 
the present, for (ii) the evolution of the BH mass function and 
accretion rate function and its connection to the observational 
"downsizing" phenomenon, for (iii) the correlations between 
black hole mass, velocity dispersion, and stellar mass of host 
galaxies, and finally, for (iv) the formation and fate of the first 
quasars and the properties of their hosts. 

An important and highly encouraging first result has has 
been that the cosmological black hole mass density pbh pre- 
dicted by our high-resolution simulation reproduces the lo- 
cally measured value and its extrapolation to higher red- 
shift (z < 2.5), inferred from integrating t he X-ray lumi- 
nosity functions ( Yu & Tremaine 2002; Ma rconi et al.l 120041 : 
IShankar et al.ll2004l) ror. more directly, the b olometric quasar 
luminosity function dHopkins et alj l2007bl) . We predict a 
steep evolution of pbh as a function of redshift, with a rise 
that is more rapid than that of the star formation rate. At z > 3, 
Pbh/p* fx (1 +z)^"' whereas there is at most weak evolution 
in the ratio Pbh/p* at z < 3. Similarly, whilst the star for- 
mation rate density p* broadly tracks the BHAR density pbh 
below z ^ 2.5 — 3, their ratio evolves steeply at higher red- 
shifts, as Pbh/p* « (1 +z)^'*. The SFR density peaks earlier 
than the BHAR density and exhibits a more gradual evolution 
with redshift compared to the BHAR. Only at redshifts below 
the peak of the BHAR, the BHAR and SFR rate densities start 
tracking each other 

Our results for the evolution of the BHAR densi ty are 
broadly consistent with constraints obtained by Hopkin s^et alj 
(2007b) from a comprehensive analysis of a large sample of 
observational data sets, which allowed them to synthesize the 
evolution of the bolometric quasar luminosity density with 
redshift and to show that the luminosity density indeed peaks 
at z ~ 2 — 3, with a sharp drop towards higher redshifts. 

We have shown that the growth of the black hole mass func- 
tion shows signs of 'anti-hierarchical' behavior, or 'downsiz- 
ing'. The high mass end of the BH mass function is estab- 
lished at comparatively high redshift and then significantly 
slows down in evolution below z 2, where only the abun- 
dance of intermediate mass BHs is still growing appreciably. 
We have found more direct evidence for black hole "downsiz- 
ing" by inspecting the distribution of accretion rates, in unis of 
the critical Eddington rate. The accretion rate function shifts 
from a narrow distribution dominated by high Eddington rates 
at high redshift to a broad distribution with a small fraction 
of Eddington accretion for low redshifts (z < 3). High mass 
black holes, forming in the high density peaks at high red- 
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Fig. 16. — Black hole merger history trees shown as a function of redshift (y-axis) and position (along the .v-axis), for the most massive black hole at z = 1 in 
the BHCosmo ran. The full tree is shown in the top left panel, while the other three panels split up the tree according to accretion rate in units of Eddington, as 
indicated by the colors. The black hole mass is given by the size of the symbols, as labeled. 



shift are built up rapidly by vigorous accretion. However, the 
effects of gas depletion and AGN feedback drive a strong de- 
cline in the accretion rate of massive black holes towards late 
times, while during these later times the peak of the accretion 
activity shifts to progressively smaller mass scales. 

The BH masses of our simulated galaxies are strongly cor- 
related with the stellar velocity dispersions and stellar masses 
of their host galaxies, and the correlations agree remarkably 
well with the local Mbh — cr and Mbh — relationships, 
over a large dynamic range. We previ ously showed with sim- 
ulations of isol ated galaxy mergers ( Di Matteo et alj [20051: 
[Robertson et al.ll2006b: Hopkins et al.ll2007ah that our AGN 
feedback prescription leads to a self-regulated BH growth that 
can explain the Mbh — cr relation. It is highly reassuring that 
the same simulation model, with an unmodified feedback effi- 
ciency e f reproduces the observed Mbh — cr in full cosmologi- 
cal simulations as well. We emphasize that the free parameter 
e / sets the normalization of the obtained relationship, but the 
slope and scatter of the relation obtained from the simulations 
are not adjustable and a non-trivial consequence of the self- 
regulated BH growth. 

Our simulation also suggests a weak evolution with redshift 
of the normalization of the Mbh ^ c relation, as Mbh/c ^ 
(1 However, we find that this evolutionary trend is 

sensitive to range of masses being probed. When we focus on 
the better resolved more massive systems, there appears to be 
some mild evolution in the slope of the BH scaling relations 



at high a and for high M» . In particular, we find a trend of 
increasing M,/Mbh with redshift, in agreement with recent 
direct estim ates of the BH to host stellar mass ratio at high 
redshift by Pen g et al] (12006V This trend is accompanied by 
an increase in the cold gas fraction in the host galaxy, and a 
smaller at a fixed stellar mass. These results are consistent 
with the suggestion that the Mbh — c and Mbh — M» relations 
are projections o f a more basic correla tion, the black hole fun- 
damental plane ( Hopkins et al.l2007a !). Mbh ct^/?"^, and the 
interpretation that gas-rich systems at high redshift produce 
more dissipative mergers that lead to more concentrated BH 
hosts. 

Our findings for the BH scaling relations appear fully 
consistent with a scenario where quasar activity is driven 
by galaxy mergers, as suggested by our si mulations of iso- 
lated merge rs (Di Matteo et al. 2005; Scringel et al. 20053 
[Robertson et al. 2006b; Ho pkins et a l. 2006a." 2007a) . While 
gas is available for star formation over a large range of mass 
scales, it can only gets to central regions of galaxies in large 
amounts as a result of the inflows that accompany major merg- 
ers. At the same time, the mergers produce spheroids, which 
together with the growth-limiting quasar feedback establishes 
the Mbh ~ o" relationship. The detailed time evolution of the 
quasar activity in an individual galaxy depends on the angular 
momentum of the gas and the disk size and morphology. It 
is clear that the limited spatial resolution of our cosmological 
simulation is a severe limitation on how faithfully this can be 
tracked in our calculation, but because the final BH masses 
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Fig. 17. — Black hole merger history tree shown as a function of redshift (y-axis) and position (.v-axis) for the first high-redshift massive black hole that forms 
in the simulation. The full tree is shown in the top left panel, while the other three panels split up the tree according to accretion rate in units of Eddington, as 
indicated by the c olor s. The black hole mass is given by the size of the symbols, as labeled. The evolution of this system is qualitatively very different from the 
one shown in Fig.[T6]Note that discreteness in the outputs for black holes that remain inactive cause the gap in the rightmost branch of this tree. 



in merger remnants are robustly predicted by our simulation 
methodology even at coarse resolution this does not seriously 
affect the quantities we study here. 

Our simulated volumes are too small to directly check if a 
population of supermassive black holes as massive the z ^ 6 
Sloan quasars can form at high redshift. However, we do 
find that our simulation model produces black holes at early 
times that spend a large fraction of their time at high accre- 
tion rates close to their Eddington rates. We therefore pre- 
dict the presence of massive black holes at high redshift. Us- 
ing high-resolution simulations of multiple mergers that were 
constructed to represent a high redshift merging tree of one of 
the most mas s ive pr otocluster regions in a (1/z^'Gpc)-' vol- 
ume, iLi et alJ (l2007h have recently made the match with the 
Sloan quasar much more explicit. They were able to show 
that the gas-rich mergers expected in this rare overdensity can 
indeed grow a BH of mass ^ 10^ Mq early enough. 

We have also found that the first supermassive black hole in 
our simulation forms in relative isolation, as a result of strong 
gas inflows and merging at the intersection of large-scale fil- 
aments. However, by following this system to lower redshift 
it turned out that it did not evolve into the most massive black 
hole today. Instead, it was overtaken in growth by a black 
hole that acquired most of its mass in major mergers in a high 
density region at lower redshift. Such a change in the rela- 
tive rank of BHs as a function of time appears quite generic, 
as we explicitly checked with a simulation of larger volume. 



This clearly complicates attempts to directly link high redshift 
progenitor systems to low redshift descendants. As we have 
demonstrated, the cosmological simulation methodology we 
introduced here provides however an excellent tool for study- 
ing the evolution of the cosmic BH evolution. In the present 
study, we focused only on the radiatively efficient accretion 
mode of AGN, which is associated with quasar activity. For 
following AGN activity also in clusters of galaxies, and hence 
down to z = in large volumes, we also need to account radio 
activity, which becomes important in very massive halos at 
low redshift. In Siiacki et al. (2007) we present an extension 
of our simulation methodology that accounts for this physics, 
and we discuss results obtained with this unified model for 
AGN feedback both for clusters of galaxies and cosmological 
boxes. 

Our cosmological approach to black hole formation enables 
us to follow the fuel for black hole activity from its ultimate 
source, the early intergalactic medium, and so link the large- 
scale structure and environments of black holes directly with 
their growth. As a result, countless different avenues of re- 
search are opened up and we plan to explore many in future 
work. For example, in Colbe rg et alj (120071) we use merger 
trees constructed for every black hole in the simulation to 
carry out a systematic study of their environments and histo- 
ries and to determine the amount of black hole growth owing 
to mergers versus gas accretion; the latter should be testable 
with upcoming gravitational wave experiments. Future issues 
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which will be addressed include the nature of AGN clustering 
and its evolution and its relation to the underlying dark mat- 
ter clustering, predictions for quasar hfetimes and luminos- 
ity functions, and how feedback from AGN may manifest it- 
self on large-scales and can be detected through the Sunyaev- 
Zeldovich effect. With our new approach we can address the 
question of how the first black hole formed and grew we will 
be able to make predictions for the ionizing background ow- 
ing to the first miniquasars, an important but currently uncer- 



tain ingredient in current models of reionization. 
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